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LONG-TERM  GOALS 

The  long-term  goal  of  this  research  is  to  construct  the  Navy’s  next-generation  global  numerical 
weather  prediction  (NWP)  models  using  new  numerical  methods  specifically  designed  for  distributed- 
memory  and  vector  computers.  To  take  full  advantage  of  distributed-memory  computers,  the  global 
domains  of  our  models  are  partitioned  into  local  sub-domains,  or  elements,  which  can  then  be  solved 
independently  on  multiple  processors.  The  numerical  methods  used  on  these  sub-domains  are  local  in 
nature,  high-order  accurate,  and  highly  efficient.  Using  these  ideas  we  are  developing  global  and 
mesoscale  atmospheric  models  that  should  be  able  to  improve  upon  the  operational  models  currently 
used  by  the  Navy. 

OBJECTIVES 

The  objective  of  this  project  is  to  construct  new  high-order  local  methods  for  the  Navy’s  next- 
generation  global  and  mesoscale  NWP  models.  The  high-order  accuracy  of  these  methods  will  ensure 
that  the  new  model  yields  the  same  accuracy  as  the  current  global  spherical  harmonics  model  and 
better  accuracy  than  the  current  mesoscale  finite  difference  model.  The  objective  is  to  achieve  this 
accuracy  while  increasing  the  geometric  flexibility  to  use  any  grid  as  well  as  to  increase  the  efficiency 
of  these  models  on  large  processor-count  distributed-memory  computers. 

APPROACH 

To  meet  our  objectives  we  explore: 

1 .  spectral  element  (SE)  and  discontinuous  Galerkin  (DG)  methods, 

2.  semi-implicit  (SI)  and  semi-Lagrangian  (SL)  time-integrators,  and 

3.  various  forms  of  the  governing  equations  in  order  to  maximize  accuracy,  efficiency,  stability,  and 
conservation  properties. 

The  power  of  SE  and  DG  methods  is  that  they  are  high-order  accurate,  like  spherical  harmonics 
methods,  yet  are  completely  local  in  nature  -  meaning  that  the  equations  are  solved  independently 
within  each  individual  element  and  processor.  Furthermore,  high-order  methods  have  minimal 
dispersion  error  which  is  an  important  property  for  capturing  fine-scale  atmospheric  phenomena  (e.g., 
tropical  cyclones,  Kelvin  and  Rossby  waves).  In  addition,  semi-implicit  (SI)  and  semi-Lagrangian 
(SL)  methods  offer  vast  improvements  in  efficiency  due  to  the  longer  time  steps  that  they  permit. 
Before  developing  a  computer  model,  it  is  important  to  first  analyze  the  best  set  of  equations  for 
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representing  the  dynamics  of  the  atmosphere.  We  have  analyzed  in  detail  various  fonns  of  the 
equations. 

WORK  COMPLETED 

NSEAM  Global  Hydrostatic  Model.  The  bulk  of  the  work  this  year  on  NSEAM  was  supporting  the 
effort  of  NRL  Monterey  to  include  the  NOGAPS  physical  parameterization  into  the  NSEAM 
dynamical  core.  For  example,  one  problem  that  arose  during  this  phase  of  the  work  was  that  the 
moisture  field  exhibited  the  so-called  “spectral  rain”  near  the  boundaries  of  the  element  edges.  A 
simple  one-dimensional  analysis  revealed  the  source  and  the  easiest  remedy  was  found  to  be  hyper¬ 
viscosity  operators  similar  to  those  used  in  current  spherical  harmonics  models.  However,  this  required 
constructing  high-order  Laplacian  operators  which  are  not  so  straightforward  to  build  in  a  finite 
element  model. 

Another  part  of  this  work  involved  constructing  the  so-called  ITCZ  (inter-tropical  convergence  zone) 
grid  which  was  used  to  study  Madden- Julian  oscillations  (MJO)  within  the  aquaplanet  experiments 
used  to  test  NSEAM  with  full  moist  physics.  This  grid  uses  high  resolution  near  the  equatorial  region 
and,  it  was  conjectured,  would  result  in  a  sharper  MJO  signal.  These  results  are  discussed  in  the  NRL 
Monterey  reports  and  a  paper  (with  Maria  Flatau)  was  submitted  to  a  journal  (see  ref.  [1]).  The  PI 
(Giraldo)  worked  on  developing  the  numerical  model  along  with  the  grids  and  the  interpolation 
stencils  required  to  output  the  fields  to  the  NOGAPS  format. 

One  portion  of  the  work  regarding  the  NSEAM  dynamics  involved  the  construction  of  a 
nonhydrostatic  DG  version  of  the  NSEAM  dynamical  core.  The  collaboration  with  Dr.  Matthias 
Lauter  (of  the  Alfred-Wegener  Institute  in  Potsdam  Germany)  through  ONR  Global  VSP  resulted  in  a 
paper  on  a  DG  barotropic  model  on  the  sphere  (see  ref.  [2]).  This  collaboration  will  continue  and  it  is 
expected  that  it  will  result  in  a  DG  baroclinic  model  by  the  end  of  2008. 

Mesoscale  Nonhydrostatic  Atmospheric  Models.  The  majority  of  the  time  this  year  was  spent  on 
analyzing  various  forms  of  the  governing  equations  of  motion  as  well  as  developing  numerical  models 
based  on  all  these  equation  sets.  Specifically  we  analyzed  the  following  equations:  Set  1  is  defined  as 
follows 

dn  R  s. 

—  +  u*V;r  +  — /r(V«u)  =  0 

dt  cv 

011  /  N  , 

—  +  u»  Vu  +  cp0V  n  =  -/(kxu)-gk  +  /jV~  u 

—  +  u*V6>=  iN20 
dt 

where  the  solution  vector  is  Exner  pressure,  velocity,  and  potential  temperature.  Set  1  is  the  equation 
set  used  in  the  U.S.  Navy’s  mesoscale  model  COAMPS  and  a  similar  fonn  is  also  used  in  the  German 
meoscale  model  LM.  The  main  problem  with  this  equation  set  is  that  it  cannot  conserver  either  mass  or 
energy. 


2 


Set  2  is  defined  as  follows 


dp 

-f-  + V»U  =  0 
dt 


fU(g>U  'j  ,  x  . 

V  •  ^ - +  PI2J  =  -/(kxU)-  /?gk  +  V  •  (jupV  u) 


5U 

dt  \  p 


where  the  solution  vector  is  density,  momentum,  and  density  potential  temperature.  Set  2  is  the  fonn 
used  in  WRF.  This  form  is  very  attractive  because  it  conserves  mass  but  not  energy. 

Set  3  is  defined  as 


^  +  v.u  =  o 

dt 


5U 

dt 


U®U 

V  P 


Ph 


=  -/(kxu)-  pgk+v»  r;; 


®  +  v.A±G,U.f; 

dt  {  p  ) 


where  the  solution  vector  is  density,  momentum,  and  density  total  energy.  This  equation  set  is  not  used 
in  atmospheric  modeling  but  is  the  equation  of  choice  in  computational  fluid  dynamics  (CFD).  This  set 
is  very  attractive  because  it  is  fully  conservative  regardless  of  whether  the  flow  is  inviscid  or  viscous. 
One  question  about  this  set,  however,  is  whether  it  can  be  coupled  to  existing  physical 
parameterization  packages  which  rely  on  potential  temperature  and  not  energy.  Using  these  three 
equation  sets  we  developed  5  different  x-z  slice  mesoscale  models  in  order  to  compare  the  spectral 
element  and  discontinuous  Galerkin  methods  and  to  see  how  these  models  behaved  under  a  series  of 
test  cases  including  sharp  front  simulations  and  nonhydrostatic  flow  over  mountains.  This  work 
resulted  in  a  peer-reviewed  article  which  was  submitted  this  year  (see  ref  [3]). 

In  addition,  we  developed  hybrid  Eulerian-Lagrangian  semi-implicit  (HELSI)  time-integrators  for  one 
of  the  SE  mesoscale  models  (for  set  1);  this  work  was  carried  out  by  Air  Force  Lt.  Tom  DeLuca  during 
his  Master’s  thesis  project  (see  refs.  [4]  and  [5]).  It  is  expected  that  this  work  will  be  submitted  to  a 
journal  and  will  continue  as  the  PhD  dissertation  of  Tom  DeLuca. 

In  collaboration  with  Marco  Restelli  (an  ONR  Global  VSP  visitor  working  on  his  PhD  dissertation 
with  Prof.  Giraldo)  the  PI  (Giraldo)  was  able  to  developed  the  first  semi-implicit  time-integrator  for 
any  DG  model.  This  result  is  very  important  because,  until  this  point,  DG  models  were  not  competitive 
with  SE  models  in  terms  of  efficiency  because  no  one  had  figured  out  how  to  construct  semi-implicit 
time-integrators  for  this  method.  This  work  has  been  submitted  to  a  journal  (see  ref.  [6]). 

Another  interesting  work  involved  the  development  of  high-order  non-reflective  boundary  conditions 
for  the  mesoscale  models.  This  work  is  being  conducted  by  Air  Force  Maj.  John  Dea  for  his  PhD 
dissertation.  The  results  are  still  quite  preliminary  but  it  looks  possible  to  be  able  to  absorb  all 
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outgoing  waves  without  any  reflections  while  using  very  few  points.  This  work  has  been  submitted  to  a 
journal  (see  ref.  [7]). 

RESULTS 

Mesoscale  Nonhydrostatic  Atmospheric  Models.  In  Ref.  [3]  we  present  five  mesoscale  models  using 
equation  sets  1,  2,  and  3,  in  conjunction  with  the  SE  and  DG  methods.  For  example,  SE1  means  that 
the  SE  method  was  used  for  set  1.  Note  that  there  is  no  DG1  because  set  1  is  not  in  conservation  form 
and  thereby  cannot  be  used  with  the  DG  method.  In  Ref.  [3]  we  show  results  for  six  test  cases  ranging 
from  rising  bubbles,  density  currents,  to  mountain  wave  problems;  for  brevity  we  only  show  two 
results  in  this  report.  Figure  1  shows  the  results  for  a  non-smooth  bubble  rising  in  an  isothermal 
atmosphere.  Figure  la  (left  panel)  shows  the  potential  temperature  perturbation  contours  after  600 
seconds.  Note  that  this  problem  represents  a  very  sharp  discontinuity  (shock  wave)  and  is  handled 
extremely  well  by  all  five  models.  Figure  lb  (right  panel)  shows  the  profiles  along  1000  meter  height 
and  here  it  is  quite  evident  that  the  discontinuities  are  quite  sharp.  Another  good  result  illustrated  here 
is  that  these  simulations  were  performed  with  semi-implicit  (for  DG3)  and  semi-implicit  semi- 
Lagrangian  (for  SE1)  methods  which  allowed  for  Courant  numbers  from  100  all  the  way  to  800  to  be 
used.  This  is  in  stark  contrast  to  Courant  numbers  of  0.5  which  are  required  by  the  explicit  time- 
integrators  (such  as  the  3ld  Order  Runge-Kutta  method  used  for  SE2,  SE3,  and  DG2  models). 


a)  b) 

Figure  1:  Potential  Temperature  for  5  meter  of  a  non-smooth  bubble  (shock  wave) 
a)  contours  after  600  seconds  and  b)  profiles  along  z=1000  meters  for  all  five  models. 
DG3  uses  the  semi-implicit  method  with  Courant  numbers  of  100  while  SE1 
uses  the  HELSI  method  with  Courant  numbers  of 800. 


In  Fig.  2  we  show  the  vertical  velocity  contours  for  two  mountain  wave  problems  where  the  dashed  red 
lines  denote  the  analytic  solutions  while  the  solid  black  lines  are  the  numerical  solutions.  Figure  2a 
(left  panel)  shows  the  results  for  the  linear  hydrostatic  mountain  wave  and  Fig.  2b  (right  panel)  shows 
the  results  for  the  linear  nonhydrostatic  mountain  wave.  Note  how  accurately  the  numerical  results 
agree  with  the  analytic  solutions. 
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Figure  2:  Vertical  velocity  contours  showing  the  analytic  solution  (dashed  red  lines)  and  the 
numerical  solution  (solid  black  lines)  for  a)  linear  hydrostatic  mountain  wave  problem  and 
b)  linear  nonhydrostatic  mountain  wave  problem.  The  results  were  obtained  with  the 
semi-implicit  DG3  with  a  Courant  number  of  1.8  and  the  HELSI SE  with  a  Courant  number  of  9. 


In  fact,  comparing  root-mean-square  errors  we  show  that  all  five  of  our  models  yield  results  that 
compare  quite  favorably  with  all  the  operational  models  currently  in  use.  Furthermore,  our  results 
show  that  the  most  important  thing  required  to  achieve  good  results  for  this  test  case  are  the  boundary 
conditions,  followed  by  the  equation  set  used.  Specifically,  set  2  (SE2  and  DG2)  is  better  than  set  1 
(SE1)  while  set  3  (SE3  and  DG3)  is  better  than  set  2.  Although  the  results  between  the  SE  and  DG 
methods  are  very  close,  it  is  clear  that  the  DG  method  outperforms  the  SE  method  in  terms  of  accuracy 
but  at  the  price  of  requiring  more  computational  resources.  Based  on  the  results  presented  in  Ref.  [3] 
one  would  conclude  that  SE3  and  DG3  are  the  best  choices  for  the  development  of  future  models. 
However,  if  one  takes  efficiency  into  consideration,  then  perhaps  the  optimal  choice  is  to  use  SE2  in 
conjunction  with  the  hybrid  Eulerian-Lagrangian  Semi-Implicit  (HELSI)  time-integrators.  Using  the 
HELSI  approach,  we  have  been  able  to  use  Courant  numbers  as  large  as  9  for  this  problem. 

IMPACT 

NOGAPS  and  COAMPS  are  run  operationally  by  FNMOC  and  is  the  heart  of  the  Navy’s  operational 
support  to  nearly  all  DOD  users  worldwide.  This  work  targets  the  next-generation  of  these  systems  for 
massively  parallel  computer  architectures.  NSEAM  and  its  mesoscale  cousins  have  been  designed 
specifically  for  these  types  of  computer  architectures  while  yielding  the  same  high-order  accuracy  as 
the  current  systems.  Additionally,  the  new  models  are  expected  to  conserve  all  quantities  such  as  mass 
and  energy  which  will  greatly  improve  the  capabilities  of  the  Navy’s  forecast  systems. 

TRANSITIONS 

Improved  algorithms  for  model  processes  will  be  transitioned  to  6.4  (PE  0603207N)  as  they  are  ready, 
and  will  ultimately  be  transitioned  to  FNMOC  with  future  NOGAPS  upgrades. 
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RELATED  PROJECTS 


Some  of  the  technology  developed  for  this  project  will  be  used  immediately  to  improve  the  current 
spectral  transform  formulation  of  NOGAPS  in  other  NRL  projects.  In  addition,  the  work  on  the 
mesoscale  models  will  help  improve  COAMPS. 
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